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ABSTRACT 

We have revisited the extended excursion set theory in modified gravity models, taking 
the chameleon model as an example. Instead of specifying their Lagrangian size, here 
we define the environments by the Eulerian size, chosen to be of the same order of 
the Compton length of the scalar field by physical arguments. We find that the Eule- 
rian and Lagrangian environments have very different environmental density contrast 
probability distributions, the former more likely to have high matter density, which in 
turn suppressing the effect of the fifth force in matter clustering and halo formation. 
The use of Eulerian environments also evades the unphysical restriction of having an 
upper mass limit in the case of Lagrangian environments. Two methods of computing 
the unconditional mass functions, numerical integration and Monte Carlo simulation, 
are discussed and found to give consistent predictions. 
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1 INTRODUCTION 

The excursion set theory (Bond et al. 1991) is a concise yet 
very successful approach (Zentner 2007) to study the non- 
linear structure formation in the standard cold dark matter 
(CDM) scenario. Starting from an (usually Gaussian) ini- 
■ tial distribution of the density perturbation field at early 
times, with an evolution model such as the spherical col- 
lapse model, 1 to map an initial density perturbation to a 
nonlinear structure (dark matter halo) at some late time, 
it can predict statistically what fraction of matter has been 
assembled into the halo at that time. It makes the physics of 
structure formation clear and simple, and was a major tool 
for studying the large-scale structure formation when large 
cosmological simulations were still beyond the capability of 
supercomputers . 

With the coming of the era of precision cosmology 
and progresses in supercomputing, iV-body simulations 
(Bertschinger 1998) have become more common nowadays 
for its ability to capture non-linear evolution without any as- 
sumptions in the evolution model and hence more accurate 
predictions than the excursion set theory. Despite this, the 
analytical results of the latter still provide valuable infor- 
mation. For example, fitting formulae of the mass function 

* E-mail: baojiu.li@durham.ac.uk 
f E-mail: tszyan.lam@ipmu.jp 

1 In this work we only consider the spherical collapse model and 
focus on the effect of different definitions of environment on the 
halo mass function in modified gravity models 



can be obtained (Sheth & Tormen 2002), the parameters 
of which are then calibrated by the numerical simulations 
(Sheth & Tormen 1999; Jenkins et al. 2001; Tinker et al. 
2008). 

In many aspects, the study of the cause of the acceler- 
ated cosmological expansion (Copeland et al. 2006) follows 
the same history of the enquiries of CDM: people first ask 
about their implications in the large, linear structures, and 
then gradually shift towards the smaller, nonlinear scales, 
in which process ever advanced techniques and tools are de- 
veloped. After an initial burst of theoretical or phenomeno- 
logical models following the first observational evidences of 
the cosmic acceleration (Riess et al. 1998; Perlmutter et al. 
1999), people have spent years trying to understand the be- 
haviour of all these models on large scales, paving the road 
that leads to a better and full understanding of the theories 
and preparing for the confrontations with future data. 

One important class of such models involves modifica- 
tions to Einstein's General Relativity (GR) on large scales 
(Jain & Khoury 2010; Clifton et al. 2011). Although these 
modified gravity models are mainly designed to explain the 
observations, many of the ideas are motivated from stud- 
ies in fundamental physics, making them very appealing. 
Clearly, because GR has been tested rigorously in the lab- 
oratories and solar system (Will 2006), any modifications 
to it must be strongly suppressed in the local environments 
and every successful modified gravity model must have some 
mechanisms to achieve this suppression to pass the first 
test. In this work we will focus on a class of models where 
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GR is modified by an additional dynamical scalar degree 
of freedom (a scalar field) which mediates a fifth force of 
gravitational strength between matter particles; the equa- 
tion of motion of this scalar makes it extremely heavy and 
therefore hard to propagate, or extremely weakly coupled 
to matter, in regions of high matter density. In both cases 
the fifth force is suppressed and therefore much weaker than 
gravity. The chameleon model of (Khoury & Weltman 2004; 
Mota & Shaw 2007) is a representative example. 

Chameleon models can have very rich phenomenology. 
In many cases, the background evolution can be indistin- 
guishable from that of the standard ACDM paradigm (see, 
e.g., Hu & Sawicki 2007; Li & Zhao 2009; Brax et al. 2008, 
for some examples). The linear perturbations on very large 
scales are unaffected by the fifth force either, because even in 
low density environments the range of the fifth force is only 
of order Mpc (Li & Barrow 2007; Li & Zhao 2009) - which 
means that observables such as the cosmic microwave back- 
ground (CMB) spectrum and integrated Sachs Wolfe (ISW) 
effect are the same as the ACDM predictions. Finally, be- 
cause of the strong suppression in the Solar system, there 
is no detectable deviation from GR locally. Consequently, 
the only place where we could expect to see effects of modi- 
fied gravity would be the nonlinear structures such as voids, 
clusters and galaxies, which are exactly the regime where 
numerical simulations are needed to make accurate predic- 
tions 2 . 

Although a number of numerical simulations have al- 
ready been done for such modified gravity models (Oyaizu 
2008; Oyaizu et al. 2008; Schmidt et al. 2009; Li & Zhao 
2009, 2010; Li & Barrow 2011; Zhao, Li & Koyama 2011; 
Li & Hu 2011; Li, Zhao, Teyssier & Koyama 2012), they are 
still at a very early stage. The reason is that the equation 
governing the scalar field is generally very nonlinear and it 
usually takes much longer to solve it than the standard Pois- 
son equation of Newtonian gravity. Performing A-body sim- 
ulations for modified gravity with very large box sizes and 
high mass/force resolutions are still a technical challenge, 
and this fact brings us back to the analytical methods, such 
as the excursion set theory. 

Unfortunately, even the application of the excursion 
set theory becomes very nontrivial in the modified gravity 
theories. The nonlinear equation essentially makes the be- 
haviour of the scalar field (and the fifth force) sensitively 
dependent on the surrounding environment, and to deter- 
mine the mapping from an initial overdensity to a late- 
time nonlinear structure we need to have knowledge of the 
environments, which have very different densities and can 
at best be described by some probability distributions (see 
Brax, Rosenfeld & Steer (2008) for an earlier study of such 
mapping). To solve this problem, Li & Efstathiou (2012) 
proposed an extension to the standard excursion set theory 
by solving the above mapping in some specific environments 
then averaging over the probability distribution of the en- 
vironments. An alternative way of viewing this approach is 
that the critical density for halo formation follows a dis- 
tribution which depends on the environment density con- 



2 Recently, Brax & Valageas (2012) has attempted to use pertur- 
bation theory to study the structure formation of modified gravity 
theories into the nonlinear regime. 



trast - in the language of excursion set theory the barrier is 
stochastic even though the collapse is described by the deter- 
ministic spherical collapse model. This stochasticity of bar- 
rier is in contrast to the one discussed in Maggiore & Riotto 
(2010); Corasaniti & Achitouv (2011) in the ellipsoidal col- 
lapse model and we refer the readers to the Appendix C 
of Paranjape, Lam & Sheth (2012a) for the discussion re- 
garding the scatter of critical density observed in numerical 
simulations and the stochasticity of barrier in the excursion 
set formalism. 

The crucial component of the extended excursion set 
model is the environment: how do we define it? Obviously, 
different definitions may give rise to different environmental 
probability distributions, and there must be some physical 
arguments to motivate the definition. As a first example to 
illustrate the idea, Li & Efstathiou (2012) define the envi- 
ronment by fixing a Lagrangian (or initial comoving) size, 
which is the simplest possibility, surrounding each proto- 
halo. As we shall discuss below, this definition has several 
drawbacks which can be cured by defining the environments 
as an Eulerian (physical) size. 

This paper is organised as follows: in § 2 we briefly re- 
view the theoretical model to be considered and summarise 
its main ingredients. We then discuss our new definition of 
environment and its motivation in § 3. § 4 is the main part 
of this work, which gives the numerical results of the uncon- 
ditional mass function using both Lagrangian and Eulerian 
definitions of environments and discusses their difference; it 
also compares the two different methods of calculation - nu- 
merical integration and Monte Carlo simulation - and finds 
good agreement. The summary and conclusions of this paper 
can be found in S 5. 



2 THE CHAMELEON THEORY 

This section lays down the theoretical framework for inves- 
tigating the effects of coupled scalar field(s) in cosmology. 
We shall present the relevant general field equations in § 2.1, 
and then specify the models analysed in this paper in § 2.2. 



2.1 Cosmology with a coupled scalar field 

The equations presented in this sub-section can be found 
in (Li & Zhao 2009, 2010; Li & Barrow 2011), and are pre- 
sented here only to make this work self-contained. 
We start from a Lagrangian density 



K 



+ V( l p)-C( V >)Cj 3 M+C s , 



(1) 



in which R is the Ricci scalar, k = 8irG with G being the 
gravitational constant, Com and Cs are respectively the 
Lagrangian densities for dark matter and standard model 
fields, <p is the scalar field and V(ip) its potential; the cou- 
pling function C(<p) characterises the coupling between ip 
and dark matter. Given the functional forms for V(ifi) and 
C(ip) a coupled scalar field model is then fully specified. 

Varying the total action with respect to the metric g a b, 
we obtain the following expression for the total energy mo- 
mentum tensor in this model: 



T ab = Vatp^blfi — gab 



+ C(<p)TT+T! b , 
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where T^ M and T^ b are the energy momentum tensors for 
(uncoupled) dark matter and standard model fields. The 
existence of the scalar field and its coupling change the form 
of the energy momentum tensor leading to potential changes 
in the background cosmology and structure formation. 

The coupling to a scalar field produces a direct interac- 
tion (fifth force) between dark matter particles due to the 
exchange of scalar quanta. This is best illustrated by the 
geodesic equation for dark matter particles 



C{ip) 



V<p, 



(3) 



where r is the position vector, t the (physical) time, <j> the 
Newtonian potential and V is the spatial derivative. C v = 
dC/dip. The second term in the right hand side is the fifth 
force and only exists for coupled matter species (dark matter 
in our model). The fifth force also changes the clustering 
properties of the dark matter. 

To solve the above two equations we need to know both 
the time evolution and the spatial distribution of ip, i.e. 
we need the solutions to the scalar field equation of motion 
(EOM) 



dVM dc(v>) _ 

V V a f H : h Pom — ; = U, 

dip dip 

or equivalently 

a dV c ff(<p) 

V V a <p H , = 0, 

d<p 

where we have defined 
Veff(<p) = V(tp) + pDMC(ip). 



(4) 



(5) 



(6) 



The background evolution of p can be solved easily given 



the present day value of pom since pdm oc 



We can 



then divide ip into two parts, ip = p> + 5tp, where (p is the 
background value and Sip is its (not necessarily small nor 
linear) perturbation, and subtract the background part of 
the scalar field equation of motion from the full equation 
to obtain the equation of motion for Sip. In the quasi-static 
limit in which we can neglect time derivatives of Sip as com- 
pared with its spatial derivatives (which turns out to be a 
good approximation on galactic and cluster scales), we find 



V p = — pdm 



dC(<p) _ dV(ip) 

-PDM + 



dV(<p) 



(7) 



dip dp dip dip 

where puM is the background dark matter density. 

The computation of the scalar field ip using the above 
equation then completes the computation of the source term 
for the Poisson equation 



[ptot + 3ptot] 

[C(<p)pDM + p B - 2V(ip)} 



(8) 



where pe is the baryon density (we have neglected the ki- 
netic energy of the scalar field because it is always very small 
for the model studied here). 



2.2 Specification of model 

As mentioned above, to fully fix a model we need to specify 
the functional forms of V(ip) and C(ip). Here we will use the 



models investigated by Li . 
with 

C{ip) = exp(7^<p), 



Zhao (2009, 2010); Li (2011), 



(9) 



and 



V{ip) = 



A 



[l-exp(-^«p)] c 



(10) 



In the above A is a parameter of mass dimension four and 
is of order the present dark energy density (ip plays the role 
of dark energy in the models) . 7, a are dimensionless pa- 
rameters controlling the strength of the coupling and the 
steepness of the potentials respectively. 

We shall choose a < 1 and 7 > as in Li & Zhao 
(2009, 2010), ensuring that V e // has a global minimum 
close to p = and d 2 V e ff(p)/dip 2 = m 2 , at this mini- 
mum is very large in high density regions. There are two 
consequences of these choices of model parameters: (1) ip 
is trapped close to zero throughout cosmic history so that 
V(ip) ~ A behaves as a cosmological constant; (2) the fifth 
force is strongly suppressed in high density regions where ip 
acquires a large mass, m 2 , 3> H 2 (H is the Hubble expan- 
sion rate), and thus the fifth force cannot propagate far. The 
suppression of the fifth force is even stronger at early times, 
and thus its influence on structure formation occurs mainly 
at late times. The environment-dependent behaviour of the 
scalar field was first investigated by Khoury & Weltman 
(2004); Mota & Shaw (2007), and is often referred to as the 
'chameleon effect'. 



3 DISCUSSION ON ENVIRONMENT 

The extended excursion set approach for chameleon models 
proposed in Li & Efstathiou (2012) differs from the original 
excursion set approach in the introduction of environment, 
which is important in the chameleon models for two reasons: 
first, the environmental density determines the critical den- 
sity for the spherical collapse inside it; second, an arbitrary 
spherical overdense region does not reside in an environment 
with any density equally likely. As a result, in the language 
of excursion set theory, the calculation of the first-crossing 
probability of the critical density curve must now be done 
in different environments and then integrated over the dis- 
tribution of the environment. 

It is therefore evident that the specification of the envi- 
ronment is crucial in the extended excursion set approach. 
In Li & Efstathiou (2012), the environment is defined as the 
follows: 

(i) it is a spherical region with a common centre as the 
considered spherical overdensity (i.e., the halo-to-be); 

(ii) it is much bigger than the halo-to-be; 

(iii) it is not too big because otherwise it will not give a 
faithful representation of the environmental density. 

As a first approximation, Li & Efstathiou (2012) defines the 
environment to have a Lagrangian radius of £ = 8/1 Mpc 
and calls this a fixed-scale environment approximation. This 
simplifies the numerical calculation and eliminates the need 
of Monte Carlo simulations. 

However, the use of the Lagrangian radius £ has certain 
drawbacks. As a first example, if the environment density is 
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high (e.g., close to the critical density for collapse), then at 
late times the sizes of the environment and the halo-to-be 
could be roughly the same, violating the above requirement 
that the environment should be much bigger than the halo- 
to-be. As a second example, if the environmental density is 
low, then it expands faster than the cosmic expansion and at 
late times can become very large in size, no longer providing 
a faithful representation of the environment. 

In Li, Zhao & Koyama (2012) it has been shown, using 
numerical simulation results, that the analytical formula for 
the fifth force used in Li & Efstathiou (2012) is quite accu- 
rate, assuming that the environment has an Eulerian (rather 
than Lagrangian) radius Q = 5 ~ 8/i _1 Mpc; see Fig. 2 there. 
As a result, a more physically reasonable definition of the en- 
vironment is its Eulerian radius. As the characteristic length 
scale in the chameleon models is the Compton length Ac of 
the scalar field (which is a function of time), a physical choice 
of the Eulerian radius £ of the environment would then be 
£ ~ O(Ac) because matter field within Ac is expected to 
affect the scalar field value at a point (this is the meaning 
of 'environment'). 

Of course, Ac evolves (increases) in time, while ( could 
either increase or decrease depending on whether the envi- 
ronmental density is lower or higher than the cosmic average. 
One can certainly choose £ ~ Ac throughout the evolution, 
but this will make the computation rather complicated. An 
alternative is to have ( ~ Ac at late times, say z — Zf where 
Zf is the formation redshift of halos. This would mean that 
at early times £ can be different from Ac , but there are two 
reasons to expect the overall effect to be small: 

(i) the fifth force only becomes important at very late 
times in most chameleon models, and for redshift z > 2 it is 
mostly negligible. At low redshifts the difference between £ 
and Ac is small; 

(ii) even at early times, the difference between £ and Ac is 
only big for overdense environments which contract during 
the cosmic history, but for these environments the fifth force 
is weak throughout the whole cosmic evolution. 

In this work we shall choose £ to be smaller than the La- 
grangian environment f; but of the same order as Ac ■ Because 
the environment is bigger than the range of the fifth force, 
its evolution could be approximated by the ACDM model. 
The Eulerian overdensity at time t can be related to the 
linearly extrapolated density contrast S(t) by the spherical 
collapse model (Bernardeau 1994; Sheth 1998): 



volume V at t: 



A NL (t) 



pV 



Ssc 



(11) 



where <5 SC « 1.676 is the critical density for ACDM at Zf = 0, 
V is the Eulerian volume, m the mass within this Eulerian 
volume and p is the cosmic background density of matter. 

The distribution of the environment density at a 
given Eulerian scale can be computed by the excur- 
sion set approach (Sheth 1998; Lam & Sheth 2008a; 
Paranjape, Lam & Sheth 2012b) using the spherical collapse 
model. Effectively Eq. (11) defines a curve B(m) in the 5- 
S plane (see below for definition of S). The first crossing of 
this environment barrier gives the value of the linear extrap- 
olated density contrast 8(t) that a spherical region contain- 
ing mass m must have in order to evolve into an Eulerian 



B(m) = 5 B , 



1 - 



-i/« s , 



(12) 



where the enclosed mass is a function of S (ses Eq. (14)). 
The mapping between m and S depends on the linear matter 
power spectrum as well as the smoothing kernel. If a power- 
law matter power spectrum P(k) with the power index n s 
is specified, then the above equation can be rewritten as 



B(S) = 5s, 



where 



c 



8/1" 1 Mpc 



3/S„. 



S 



(3 + n 3 )5 s 



(13) 



(14) 



S(m) = S(0 = k 2 P(k)W 2 (k0dk, 

with W(k£) the filter function and £ the Lagrangian radius 
of the filter so that m(£) = |7r£ 3 p. erg is given by S = a 2 . 
with £ = 8h~ Mpc. In particular for n 3 — —1.2 this barrier 
becomes linear in s (approximating Ssc = 5/3) and make 
the analytical analysis easier. In what follows we shall use 
( — 5/i _1 Mpc as suggested by Li, Zhao & Koyama (2012). 
One can make a similar analogy for Lagrangian environ- 
ment, only in this case the barrier is vertical (see figure 1) 
and there is one and only one crossing of this Lagrangian 
barrier. 

In the language of excursion set theory, then, the first 
crossing probability of the moving barrier -6(5*) in [5*, S+dS'], 
denoted by P env (S)dS , is the probability that an arbi- 
trary point (the centre of a halo-to-be) is in an environ- 
ment whose linearly extrapolated density contrast falls into 
[B(S), B(S + dS)] and this environment will evolve into a 
region having an Eulerian radius £ at zf, the total mass en- 
closed in this environment is m(S + dS) ^ M env ^ 771(5*). 
Both the environmental density and environmental mass are 
important for the discussions below. 



4 UNCONDITIONAL MASS FUNCTIONS 

As a first application of the idea described above, in the rest 
of the paper we shall study the (unconditional) mass func- 
tion of the chameleon models, using both numerical integra- 
tion and Monte Carlo simulations, assuming uncorrelated 
steps in the framework of the excursion set formalism. The 
analysis of correlated steps and the associated conditional 
mass function and halo bias will be discussed in Lam & Li 
(2012). 



4.1 Eulerian versus Lagrangian environments 

Fig. 1 illustratively demonstrates the difference between 
the two definitions of the environment (Eulerian and La- 
grangian). Here, the dotted horizontal line denotes the crit- 
ical density S sc ~ 1.676 for a ACDM collapse, and the mid- 
dle one of the dashed vertical lines represents the value of S 
corresponding to a Lagrangian radius £ = 8/i~ 1 Mpc, which 
was used in Li & Efstathiou (2012) as the definition of the 
(Lagrangian) environment. Consider now two random walk 
trajectories (the blue and red curves): they cross the middle 
vertical dashed line at the two dark grey filled triangles, but 
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Figure 1. (Colour Online) An illustration of the Eulerian versus 
Lagrangian environments. See the text for a detailed description. 
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Figure 2. The probability distribution, P (tWv), of the Eulerian 
(solid curve) and Lagrangian (dashed curve) environments, nor- 
malised to unity. 



do not cross the dotted horizontal line before that. Those 
two triangles correspond to two very different values of 8 env 
(1.3 and —0.4) and thus the two random walks represent 
structure formations in very different (Lagrangian) environ- 
ments. 

The Eulerian environment, on the other hand, is de- 
fined such that its Eulerian radius today is 5/i _1 Mpc and it 
is represented by the solid curve. The random walks cross 
this curve at the two grey filled circles, corresponding to 
5 cn v ~ 1-25 and 0.25 respectively: the environments are less 
different if they are Eulerian! In particular, the structure 
formation relevant for the red random walk now takes place 
in an overdense rather than underdense environment. 

Note that the two filled circles correspond to Lagrangian 
radii of £ = 10.3 and 5.43/i~ Mpc (the left and right verti- 



cal dashed lines) respectively, while the two filled triangles 
correspond to Eulerian radii of f = 3.65 (dash-dotted curve) 
and 8.9/i _1 Mpc (dash-dot-dot curve) respectively. It is clear 
that the Eulerian environment of the blue random walk con- 
tains about seven times more matter than that of the red 
random walk, but if one uses Lagrangian definition then 
both environments contain the same amount of matter. 

The above result has implied that the Eulerian and La- 
grangian environments should have different probability dis- 
tribution functions (PDFs). We then need to compute the 
PDF of <5 C nv, P(5 env ), for these environment definitions: the 
Lagrangian density contrast follows the Gaussian distribu- 
tion and the corresponding P(8 cnv ) is given by Eq. (40) of 
Li & Efstathiou (2012), while the Eulerian density can be 
considered as a mixture of Lagrangian density from different 
smoothing scales: a dense Eulerian environment is evolved 
from a bigger Lagrangian batch but the opposite is true for 
underdense Eulerian environments. 

There are different methods to calculate P(S onv ) for the 
Eulerian environment: one can either obtain it by solving 
the first crossing distribution of the barrier in Eq. (12) (see 
Zhang & Hui (2006); Lam & Sheth (2009) for uncorrelated 
steps; Maggiore & Riotto (2010); Corasaniti & Achitouv 

(2011) ; Paranjape, Lam & Sheth (2012a); Musso & Sheth 

(2012) for correlated steps), or using analytical expressions 
such as the log-normal distribution or the expressions given 
in Lam & Sheth (2008a, b) for the evolved nonlinear den- 
sity contrast and combining that with Eq. (11) to obtain 
P(8 C nv) (recall that S cnv is the linearly extrapolated den- 
sity contrast at the Eulerian environment). In particular, 
for power-law matter power spectrum the barrier is given by 
Eq. (13) and applying the first crossing probability approx- 
imation in Lam & Sheth (2008a) (see Lam & Sheth (2009) 
for an explanation for this approximation), the distribution 
of S cnv is 



P(<5cnv) = 



,/2 



x exp 



1 + (w - 1) 

/r 



5c 



„/2-l 



2 (i - s eav /8 c y 



(15) 



where f3 = (£/8) ' c /a 8 , ui — 5 c ^y, and 7 is the logarithmic 
derivative of the density fluctuation variance w.r.t. m: 

dlnS 1 7i s +3 , . 

dmm 3 

The case where n s — —1.2 is of special interest, since not 
only the barrier B(S) is linear in S (if one set 8 C = 5/3), the 
first crossing probability approximation above is also exact 
(see, for example, Lam & Sheth 2009, for derivation). The 
associated distribution of 5 en v becomes 

£1/2 



P{5 eav ) = 



1 



,/3/2(l-« env /i c ) 



(17) 



/27 (1 - <W5 C ) 3 / 2 
The two environment density distributions for the case 
of power-law power spectrum with n s = —1.5 are shown 
in Fig. 2, which shows the PDF for the Eulerian environ- 
ment peaks at bigger 5 cnv than that for the Lagrangian 
environment. The results were computed numerically using 
the method described in Zhang & Hui (2006) and used in 
Li & Efstathiou (2012). Note that 

(i) If S cnv is close to S sc the PDF is smaller for the Eule- 
rian environment, because for such an environment to evolve 
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Figure 3. The critical density for the spherical collapse at zj = 0, 
as a function of the linearly extrapolated environmental density 
contrast <5 en v and S. 



force is strongly suppressed so that the collapse is governed 
by Newtonian gravity only and S c ~ S S c', 

(ii) as 8 e nv decreases, the fifth force becomes stronger in 
general and enhances the matter clustering. This means that 
the required critical density to collapse at Zf is lower; 

(iii) very small values of S correspond to the very big 
smoothing scales. Such regions have efficient self-screening 
of the fifth force irrespective of their environment, and hence 
5 C (S, Senv) ~ S sc for arbitrary <5 cn v for small S; 

(iv) large values of S correspond to small smoothing 
scales, which does not have strong self-screening of the fifth 
force, and hence 6 C (S, 8 env ) depends sensitively on <5 en v 

In practice, once we know <5 cnv , we can fully determine 
the collapse criteria for halos of arbitrary size in this en- 
vironment, and therefore the first crossing probability den- 
sity, f(S\S C nv), across the barrier S C (S, <5 en v). Here | denotes 
it is the conditional first crossing probability that the ran- 
dom walks having first crossed B(S) given in Eq. (13) at 
B(S) — S cnv . One then only needs to average over the 
first crossing probability distribution of <5 on v (which we have 
found in the previous subsection) to find the averaged first 
crossing probability f(S). 



into an Eulerian radius of 5h 1 Mpc today, it must have 
an extremely large Lagrangian size (£ > 8ft _1 Mpc), which 
is a very rare event; Note that the other reason is due to 
the approximation formula Eq. (11) - it maps 5u n — > S c to 
<5 n i — > oo. Hence p(5 cnv ) must go to zero for S cnv ^ S c . 

(ii) The PDF for Lagrangian environment peaks at S cnv ~ 
0, because at any value of S (say S = 0.64 for £ = 8/i _1 Mpc, 
middle vertical dashed line of Fig. 1) the value of the random 
walk position satisfies a Gaussian distribution; 

(iii) The PDF for Eulerian environment peaks at S cnv ~ 
0.8 > 0, because the random walk is crossing a decreasing 
barrier. This means that if one integrates over the PDF of 
environments then more contributions come from higher- 
density environments under the Eulerian definition; 

(iv) The PDF is lower for small 5 cnv (< —1.5) under the 
Lagrangian definition, because underdensities with linearly 
extrapolated density contrast 8 < —1.5 and Lagrangian size 
£ = 8/i _1 Mpc will have evolved into very large Eulerian sizes 
today (£ 3> 5/i _1 Mpc), and this is rare event; in contrast, 
to have an Eulerian size of 5/i _1 Mpc today such undensities 
can be quite small in initial sizes whose r.m.s fluctuation is 
bigger and hence the probability is not negligible. 



4.2 Spherical collapse in a given environment 

The spherical collapse history of an initial overdensity in a 
given environment of linearly extrapolated density contrast 
5 e nv has been discussed in detail in Li & Efstathiou (2012), 
to which interested readers are referred, and here we will 
only give a qualitative description of the results. 

Fig. 3 shows the critical density S C (S, 8 cnv ) for an initial 
overdensity to collapse at Zf = 0, as a function of both S and 
<5env As the fifth force is always attractive, 8 C (S, 5 cnv ) < S sc 
where S sc is the critical density if the fifth force vanishes. 
Furthermore, 



4.3 Averaging over environmental distribution 

As in Li & Efstathiou (2012), the final first crossing proba- 
bility, which is related to the unconditional mass function, 
is calculated by making the environmental average: 



f(S) 



env 5 



(18) 



where the upper limit of the integral is S sc because S cnv ^ S sc 
(c.f. Fig. 2) because by definition an environment has not 
collapsed to form a halo by zj. In the case of Lagrangian en- 
vironment, P(S cnv ) is defined so that it is identically zero for 
Scnv ^ <5 S c, as shown in Eq. (40) of Li & Efstathiou (2012). 
For the Eulerian environment proposed in this work, the Eu- 
lerian barrier always lies below <5 SC for all S > hence this 
upper limit is valid by construction. 

The mass function, dn/dM, is related to f(S) by 



dn 
AM 



AM 



dS 



AM 



AM, 



(19) 



where n is the number density of halos and M the halo mass. 
In the literature, an alternative quantity which is often used 
is the collapsed mass fraction vF(y), with v = S^ c /S. vF(v) 
is the fraction of matter in collapsed objects per logarithmic 
interval of v and satisfies 



F{v)Av = f(S)AS 



vF(v) = Sf(S). 



(20) 



(i) in high-density environments (<5 e 



the fifth 



In the following subsections, we shall calculate the quan- 
tity vF(y) using two methods: the numerical integration of 
Eq. (18) and Monte Carlo simulations. We focus on the ex- 
cursion set with uncorrelated step in the current study but 
would like to point out that Eq. (18) applies in both cor- 
related and uncorrelated steps calculation - in the case of 
correlated steps one must take into account the fact that 
the condition in /(^^env) is the random walk first crosses 
the environment barrier at S anv , or in other words the ran- 
dom walk is non-Markovian. It complicates the analytical 
calculation and it will be discussed in Lam & Li (2012). In 
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Figure 4. (Colour Online) Upper left panel: the collapsed mass fraction uF{v) for the ACDM and chameleon model with a = 10 -6 , £t = 
0.5; for each model the environmental average (see the text) is performed assuming Lagrangian and Eulerian environments respectively 
(see the legends), and the solid curve is the analytic result for ACDM. Upper right panel: the same as the upper left panel, but here 
only the results for a chameleon model with a = 10~ 5 , [i = 0.5 (the symbols) are shown; the analytical results for ACDM is shown as a 
solid curve for comparison. Lower left panel: the relative difference of vF{y) between the chameleon model with a = 10 -6 and ACDM. 
Lower right panel: the relative difference of vF[y) between the chameleon model with a = 10~ 5 and ACDM. Note that the vertical axis 
is ln(10)^-F(i^) rather than vF(y). 



the present case where the random walk has uncorrelated 
steps (it is Markovian) , the conditional probability does not 
depend on the history of the walk prior to 8 cnv . 

In what follows we apply two approaches to evaluate 
the halo mass function in modified gravity models. The 
numerical integration approach is similar to the previous 
work by Li & Efstathiou (2012) which applies the method in 
Zhang & Hui (2006); the Monte Carlo simulation approach 
follows the variation of density contrast as a function of 
smoothing scale and we keep tracks of the first crossing of 
the environment barrier as well as the first crossing of the 
consequent modified halo formation barrier. 

4-3.1 Numerical integration 

As numerical examples, we have calculated vF(v) as a func- 
tion of v for three models - GR 3 and two chameleon models 
with fi = 0.5, a = 10~ 6 and n = 0.5, a — 10 -5 - all having 
the same background cosmology. A power-law matter power 
spectrum with index n s = — 1-5 is assumed. For each model 
we have done the calculation assuming Lagrangian and Eu- 
lerian environments respectively, and the results are shown 
in Fig. 4. 

Assuming spherical collapse, there is no environmental 

3 We will consider GR with a cosmological constant which drives 
the accelerating expansion of the Universe. So in what follows we 
will use the words 'GR' and 'ACDM' interchangeably. 



dependence in the ACDM result, which means that the pre- 
diction of vF(y) should be the same whether one uses the 
Lagrangian or Eulerian definition of environment. This has 
been confirmed in the upper left panel of Fig. 4, which serves 
as a check of accuracy of the numerical computation. For 
comparison we have also plotted the exact analytic solution 
for the ACDM model (solid curve) 




In Fig. 2 we have seen that the PDF for Eulerian envi- 
ronment peaks at higher <5 cn v than the PDF for Lagrangian 
environment. Because higher environmental density means 
stronger suppression and therefore weaker effect of the fifth 
force, we would expect that the deviation from ACDM pre- 
diction of vF(y~) be smaller if we use the Eulerian instead of 
Lagrangian definition of environment. This is confirmed by 
Fig. 4, which shows that the former gives smaller vFiv) for 
medium and large v (halos of medium and large sizes). 

For small values of v, on the other hand, vF(y) is bigger 
if one uses Eulerian environment, because a smaller fraction 
of matter has fallen into medium and large collapsed ha- 
los compared to the case of Lagrangian environments, and 
therefore more matter is left to form the isolated small halos. 

Note that one disadvantage of using a (fixed-sized) La- 
grangian radius, as in Li & Efstathiou (2012), is that all the 
environments have exactly the same initial size and so con- 
tain the same amount of matter. In particular, the biggest 
halo to form in any given environment cannot contain more 
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Figure 5. (Colour Online) Comparison of the mass fraction func- 
tions obtained from numerical integration above and Monte Carlo 
simulations. See the legends for the details. 
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matter than that is contained in the environments, and this 
places an upper limit of the halo mass, as is shown in Fig. 4. 
If one uses Eulerian environment this problem is solved be- 
cause the environment can be arbitrarily large in size. 

4-3.2 Monte Carlo simulations 

In this subsection we compare the first crossing distribution 
for different models with or without the chameleon fifth force 
using Monte Carlo simulations. It is customary to use the 
variance in Eq. (14) with a tophat window function to relate 
S and the smoothing scale R although strictly speaking it 
is not fully consistent since random walks with the tophat 
window function induces correlated steps (or one needs to 
use the sharp fc-space window function - which does not 
have a well defined enclosed mass - to obtain uncorrelated 
steps) (Bond et al. 1991). 

In the Monte Carlo simulations a sample of random 
walks is generated following the procedure in Bond et al. 
(1991). In hierarchical models, the variance of the density 
field is a monotonic function of smoothing scale. Hence vari- 
ance of the density field S, the smoothing scale as well as 
the total mass enclosed are interchangeable quantities. For 
the case of ACDM a constant barrier with S sc = 1.676 is 
assumed. For chameleon models the environment density 
Jenv is recorded at the corresponding scales: for Lagrangian 
environments, the environment density 5 cnv ,L is recored at 
£ = 8/i _1 Mpc; for Eulerian environment c>env,E is the height 
of the random walk where it first crosses the barrier B(S) 
in Eq. (12). The halo formation criterion S C (S, <5 en v) is deter- 
mined using this environment density <5 env (see Fig. 3). We 
then follow the random walk until it first crosses S C (S, 5 en v) 
and record the associated value of S. 

Note that an advantage of using an Eulerian scale as 
environment is that the Eulerian barrier B(S) by definition 
always lies below the constant barrier - hence it is impos- 
sible to reach 8 SC before crossing B(S) and an environment 
always contains more matter than the halo inside it does. On 



Figure 7. (Colour Online) Comparison of the effect of mass 
function for different environment definitions in the chameleon 
models. Results are obtained from Monte Carlo simulations with 
ACDM power spectrum where ratios of mass function are taken 
with respect to that of Eulerian environment ( = 5h~ 1 Mpc 
(red solid curve in figure 6). Three other chameleon models are 
considered: Lagrangian environment with £ = 8h~ 1 Mpc (green 
boxes); Eulerian environment with £ = Sh~ 1 Mpc (blue dotted 
curves) and with £ = 10h _1 Mpc (magenta dashed curves). Left 
panel: chameleon model with fj, = 0.5, a = 10 — 6 ; right panel: 
fj, = 0.5, a = 10 -5 . Only errors bars for Lagrangian environment 
with £ = 8fe _1 Mpc ( green boxes) are shown for clarity. 

the other hand, if one uses Lagrangian environments, it is 
possible, though rare, that the random walk reaches S sc at a 
Lagrangian scale bigger than the predefined Lagrangian en- 
vironment radius £: this is just the problem we discussed at 
the end of Sect. 4.3.1, and in this case we assume that the 
evolution of this huge overdensity follows the background 
environment (8 cnv — 0). For £ = 8/i _1 Mpc it corresponds to 
first crossing across 8 SC beyond \g(u) ~ 0.63. 

Fig. 5 compares the output of the Monte Carlo simula- 
tions to that from the numerical integration described in the 
previous subsection. For this purpose we have only plotted 
the ACDM model and a chameleon model of /j, = 0.5, a = 
10~ 6 , adopting the Eulerian definition of the environment 
with a radius of 5h~ Mpc. The predicted mass functions 
from the two methods agree with each other reasonably in 
both models and, in the former case, they are consistent 
with the analytical result as well. 

As Monte Carlo simulations with uncorrelated steps are 
generally much faster than the numerical integration, we 
shall use the former to test the effects of different physical 
parameters. Fig. 6 and 7 show the dependency of the col- 
lapsed mass fraction vF(y) on the definition of environment, 
the initial power spectrum and the model parameter a, from 
which we can see that 

(i) increasing the radius of the Eulerian environment re- 
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Figure 6. (Colour Online) Comparison of mass function obtained from Monte Carlo ssimulation. Left column: chameleon model with 
fj, = 0.5, a = 10 -6 ; right column: fj, = 0.5, a = 10 — 5 . The upper row shows results using ACDM power spectrum while the lower row 
shows results with power-law power spectrum. In each panel the black histogram shows the GR + ACDM mass function. Two choices of 
environment are chosen for the chameleon models: Eulcrian with f = 5/i~ 1 Mpc (red solid curves), Lagrangian with £ = 8h _1 Mpc (green 
boxes). Only error bars for the red solid curves are shown in the lower panel for clarity. 



suits in more big halos because it necessarily means that the 
PDF of the initial environment density S cnv shifts towards 
lower values of S cnv , making the formation of such halos more 
strongly affected by the fifth force. This dependence on the 
Eulerian environment radius is subdominant compared to 
the modification to the GR case; 

(ii) the effect of the Eulerian environmental radius on the 
abundance of small halos is weaker; 

(iii) switching to a Lagragian environment results in a 
distinctive drop in intermediate mass halo - this decrease 
corresponds to the Lagrangian scale chosen: beyond lg(z/) = 
0.63 our Lagrangian environment chameleon models reverts 
to the GR mass function; 

(iv) the results using a ACDM matter power spectrum are 
similar to the power-law power spectrum with n B = —1.5; 

(v) increasing a weakens the chameleon effect which sup- 
presses the fifth force in high matter-density regions, and 
thus producing more large halos. 



5 DISCUSSIONS AND CONCLUSIONS 

For modified gravity models which reduce to Newton grav- 
ity in high density environments such as the Solar system, 
while at the same time have significant deviations from GR 
on Mpc scales, the environment is often an important con- 
cept in both the theory and the techniques used to analyse 
it. In this paper we extended the work of Li & Efstathiou 
(2012) on estimating the modification of halo abundance in 
chameleon models within the framework of the excursion set 
approach by considering the possibility of defining environ- 
ments according to their Eulerian, rather than Lagrangian, 
sizes. Being the physical size, the Eulerian size changes with 
time and describes the dynamical nature of the true environ- 
ment. By choosing the Eulerian size to be of the same order 
as the Compton length of the scalar field, the single length 
scale in the theory, we have a better motivated definition of 
the environments. 
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Of course, the exact value of the Eulerian size £ is still a 
free parameter which can be tuned to match the simulation 
data. Alternatively, one can make certain approximations in 
the calculation to obtain the functional form of the mass 
function and then calibrate the parameters using simula- 
tions. Another possible solution is by applying correlated 
steps in the excursion set approach (see Lam & Li (2012)) - 
dramatic fluctuations between similar smoothing scales are 
unlikely with correlated steps and may be able to evade this 
ambiguity in the choice of Eulerian size. 

The Lagrangian definition of environments also suffers 
from a limitation, namely that all environments, having the 
same initial comoving size, contain the same amount of mat- 
ter. This effectively sets an upper limit of the halo mass 
that can be studied. The Eulerian definition of environment 
solves this problem because here in principle the environ- 
ment can be infinitely large, and halos of any mass can form 
within it. 

In our theoretical framework, the effect of different 
definitions of the environment enters into the calculation 
through the change in the PDF of <5 cnv , which describes 
the distribution of density contrast in the environment sur- 
rounding halo. Since the strength of the fifth force depends 
on Senv, modifications in p(S cnv ) alters the halo abundance 
for different environment defintions. In particular we found 
that it is more likely to have high values of 8 eU v in Eulerian 
environment. In this case the fifth force is suppressed and 
the net result is smaller deviation of the mass function from 
ACDM predictions, especially for the big and medium-sized 
halos. We have verified this using two methods: numerical 
integration and Monte Carlo simulations, and both methods 
agree with each other. 

To briefly summarise, this work emphasises the impor- 
tance of environment definition in the study of structure 
formation in modified gravity theories, and lays down the 
formalism for applying a physically motiviated Eulerian en- 
vironment in the framework of the excursion set approach. 
This framework could easily be generated to other types 
of modified gravity theories, such as the envrionmentally 
dependent dilation (Brax et al. 2010, 2011) and the sym- 
metron (Hinterbichler & Khoury 2010; Davis et al. 2011) 
modes. Within this framework we can easily analyse other 
quantities of interests, such as the halo bias, voids statistics 
and merger tree. We can also derive approximate analytical 
expression for the halo mass function and calibrate the pa- 
rameters with the numerical simulations. These will be the 
topics of future works. 
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